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Wc report rncasurernents of transverse inoineuturn pt spectra for ten event multiplicity classes of 
p-p collisions at ^/s — 200 GeV. By analyzing the multiplicity dependence wc find that the spectrum 
shape can be decomposed into a part with amplitude proportional to multiplicity and described by 
a Levy distribution on transverse mass mt, and a part with amplitude proportional to multiplicity 
squared and described by a gaussian distribution on transverse rapidity yt- The functional forms of 
the two parts are nearly independent of event multiplicity. The two parts can be identified with the 
soft and hard components of a two-component model of p-p collisions. This analysis then provides 
the first isolation of the hard component of the pt spectrum as a distribution of simple form on yt- 
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I. INTRODUCTION 

The structure of the inclusive pt spectrum from rel- 
ativistic nuclear collisions is affected by several aspects 
of collision dynamics and by the final-state hadroniza- 
tion process. Comparisons of p-p, d-Au and Au-Au pt 
spectra at RHIC suggest that a form of color-deconfined 
matter has been created in Au-Au collisions 1, 2]. Parti- 
cle production mechanisms which could determine spec- 
trum structure include soft parton scattering followed by 
longitudinal or 'string' fragmentation 3] and hard par- 
ton scattering followed by transverse fragmentation 
Other mechanisms could be significant. The structure of 
the Pt spectrum at some achievable level of precision may 
therefore be complex. A summary of efforts to unfold and 
interpret the structure of inclusive pt spectra from ISR 
to Fermilab and SPPS energies in the context of jet phe- 
nomenology and QCD (quantum chromodynamic) theory 
is provided in ^5]. 

At RHIC energies hard parton scattering is expected 
to dominate the spectrum at larger pt and to be signif- 
icantly modified in A-A collisions (jet quenching) 
But how does hard scattering contribute at smaller pt? 
How does it interact with thermal or 'soft' particle pro- 
duction? Is there an 'intermediate' pt region with its 
own unique production mechanisms? Those issues re- 
main unresolved after much theoretical speculation and 
experimental measurement and provide a context for the 
present analysis applied to high-statistics pt spectra from 
ten multiplicity classes of p-p collisions. The multiplic- 
ity dependence offers new access to underlying particle 
production mechanisms. 

Pt spectra from relativistic nuclear collisions are con- 
ventionally modeled by the power-law function a form 
suggested by measured jet systematics and perturbative 
QCD (pQCD) expectations. At larger pt the spectrum is 
expected to tend asymptotically to the power-law form 
p^" 9] . The strict power-law form is then generalized to 
the function A/{1 +pt/poy\ having the expected pQCD 
dependence at larger pt but transitioning to an approxi- 
mate Maxwell-Boltzmann form at smaller pt, consistent 
with expectations for thermal particle production. Al- 
though the power-law function has been previously ap- 
plied to p-p data with apparently good fit quality (x^ 
within expected limits) it has not been tested with the 
precision of recently-acquired RHIC p-p data. One can 
question the validity of its underlying assumptions. For 
instance, why should a single model function adequately 
describe spectra which may represent a mixture of several 
particle production mechanisms? 

Alternatively, a model function can be formulated 
in terms of the two-component model of nuclear col- 
lisions which identifies 'soft' p-p collisions with 
no hard parton scatter and 'semi-hard' collisions with 
at least one significant parton scatter (i.e., producing 



distinguishable hadron fragments). According to the 
two-component model the minimum-bias distribution on 
event multiplicity Uch can be decomposed into separate 
negative binomial distributions (NBD) identified with 
soft and semi-hard event types. We then expect the frac- 
tion of events with a hard parton collision to increase 
monotonically with selected event multiplicity Uch- Vari- 
ation of Pt spectra with Uch could then provide a basis for 
isolating soft and hard (and possibly other) components 
of inclusive spectra on a statistical basis, where the hard 
spectrum component refers to the fragment pt spectrum 
for hard-scattered partons, and the soft component is the 
Pt spectrum for 'soft' particle production. 

In this analysis we first test the ability of the conven- 
tional power-law model function to represent the data. 
We then reconsider the data with no a priori assump- 
tions. We attempt to describe all spectrum structure 
with the simplest algebraic model required by the data 
{e.g., 'simple' in terms of parameter number and func- 
tional forms - cf. Eq. and Sec. IXI|I and then to as- 
sociate the model elements with possible particle pro- 
duction mechanisms. We adopt two new analysis tech- 
niques: 1) We introduce transverse rapidity yt 0, [3 
as an alternative to pt. yt has the advantage that spec- 
trum structure associated with hard parton scattering 
and fragmentation is more uniformly represented on a 
logarithmic variable: yt corresponds to variable = 
ln(pparton/Pfragment) Conventionally used to describe par- 
ton fragmentation functions in elementary collisions |l3j| . 
A simple description of soft particle production is not 
compromised by the choice of transverse rapidity. 2) We 
introduce the running integral of the yt spectrum, which 
substantially reduces statistical fluctuations relative to 
significant structure and therefore improves the precision 
of the analysis. 

In this paper we present high-statistics pt spectra for 
ten multiplicity classes from p-p collisions at y/s = 200 
GeV. We use the conventional power-law model function 
to fit those spectra and assess the quality of that descrip- 
tion. We then construct running integrals of the spectra 
on yt and define a reference function common to all Uch 
values and based on the Levy distribution. We use that 
reference to extract difference spectra which contain the 
ricTi-dependent parts of the spectra in a more differential 
form. We find that the difference spectra have a simple 
structure: the major component is well-described by a 
gaussian distribution with fixed shape and with ampli- 
tude (relative to the reference) linearly proportional to 
the particle multiplicity. To simplify presentation we ini- 
tially describe approximate relationships and optimized 
parameters without errors. We then return to a compre- 
hensive discussion of the parameter system and its errors 
and consistency in Sec. IVII This analysis is based on p- 
p collisions at ^/s = 200 GeV observed with the STAR 
detector at the Relativistic Heavy Ion Collider (RHIC). 
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II. pt AND yt SPECTRA 

Data for this analysis in the form of inclusive pt 
spectra for unidentified charged particles were obtained 
from non-single-diffractive (NSD) p-p collisions at -^i = 
200 GeV triggered by a coincidence of two beam-beam 
counters (BBC) in 3.3 < hi < 5 Q- Charged particles 
were measured with the STAR Time Projection Chamber 
(TPC) and Central Trigger Barrel (CTB) [l3|. Particle 
momenta were determined with a 0.5 T magnetic field 
parallel to the beam (z) axis. Primary charged parti- 
cles were represented by TPC tracks falling within the 
acceptance for this analysis - 2tt azimuth, pseudorapid- 
ity I77I < 0.5, and 0.2 < Pt < 6 GeV/c - and satisfying 
track cuts described in The observed particle mul- 
tiplicity in the acceptance is denoted by hch, whereas 
the corrected and pt-extrapolated true event multiplicity 
is denoted by rich- From 3 x 10^ NSD events individ- 
ual Pt distributions were formed for 10 primary-particle 
multiplicity classes indexed by the observed multiplicity: 
hah G [I,-- - ,8,9+10,11 + 12]. 

To eliminate backgrounds from event pileup each TPC 
primary-track candidate was independently required to 
match a CTB/trigger timing requirement (100 ns, match- 
ing efficiency 94%, false-coincidence background 2%) and 
project to the beam line within 1 cm transverse distance 
of closest approach. No other vertex requirement was 
applied to the primary tracks. The event-vertex z posi- 
tion was estimated by the arithmetic mean z of projected 
track z for all CTB-matched primary tracks in an event. 
Events with |z| < 75 cm were accepted for further analy- 
sis. The event vertex was not included in primary-track 
Pt fits. That procedure eliminated pileup-event tracks, 
selected those events well-positioned relative to the TPC 
and minimized correlations of individual track pt and pt 
spectrum shape with event multiplicity or event trigger- 
ing not related to collision dynamics. 

The resulting pt spectra were corrected for track- 
ing efficiency, backgrounds and momentum resolution. 
Tracking acceptance and efficiency on {pt^rj, z) and back- 
grounds were determined by embedding Hijing events in 
data events with at least one empty bunch (so-called 
abort-gap events). The same fractional correction was 
applied to all multiplicity classes. The correction fac- 
tor was 1.45 at 0.2 GeV/c, falling to 1.2 at 0.5 GeV/c 
and thereafter smoothly to 1 at 6 GeV/c. Efficiency- 
and acceptance-corrected (but not pt-extrapolated) spec- 
tra integrate to multiplicity n^^ = (1.35 ± 0.015) rich, 
while the corrected and pt-extrapolated per-event spec- 
tra integrate to 'true' multiplicity rich = (2.0 + 0.02) hch- 
The errors reflect the spectrum-to-spectrum relative nor- 
malization uncertainties most relevant to this differential 
analysis. The normalization uncertainy common to all 
spectra is about 10%. 

In Fig. n (left panel) corrected and normalized per- 
event Pt spectra are plotted as points in the form 
^/n-ch^/ptdn/dpt for ten multiplicity classes, offset by 
successive factors 40 (except for hch = 1 at bottom). 
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FIG. 1: Corrected and normalized charged-particle spec- 
tra on transverse momentum pt (left) and transverse rapid- 
ity yt (right) for 10 event multiplicity classes, displaced up- 
ward by successive factors 40 relative to hch = 1 at bottom. 
Solid curves represent reference function Us/rich ■ So{yt) {cf. 
Sec Ai-V a . Dotted curves are spline fits to guide the eye. 

Parentheses for ratio prefactors of spectrum densities in 
the form dn/dx are omitted to lighten notation. In other 
cases ratio prefactors are separated from densities by a 
dot. Corrected and extrapolated spectra normalized by 
rich all integrate to unity in the sense of Eq. Q for 
Pt or yt, with the integration limit ^00. In Fig. ^ 
(right panel) equivalent spectra on transverse rapidity 
are plotted. Hard parton scattering leading to trans- 
verse fragmentation may be better described on trans- 
verse rapidity yt = In {{mt + pt) / ttiq} , with transverse 
mass mt = ^ p\ -\- rriQ and pion mass m^r assumed for 
mo: yt — 2 ^ pt ^ 0.5 GeV/c and yt = 4.5 ^ ~ 6 
GeV/c. The solid curves iig/rich • Sq provide a visual ref- 
erence for the data, risihch) and Sa{pt or yt) are defined 
below, and function is by definition independent of 
hch- 

III. POWER-LAW ANALYSIS 

The power-law function is the conventional model func- 
tion applied to pt spectra from rclativistic nuclear col- 
lisions 1^]. Said to be 'QCD-inspired,' the function 
A/{1 + pt/po)^ goes asymptotically to p^" at large pt 
(hence 'power-law') and approximates an exponential at 
small Pt- The argument supporting the power-law func- 
tion assumes that pt spectra at larger collision energies 
can be modeled with a single functional form. In this 
part of the analysis we test that assumption. The pt 
spectra for ten multiplicity classes in Fig. ^ were fitted 
with the three-parameter power-law model function de- 
fined above. Parameters A, po and n were independently 
varied to minimize for each multiplicity class (in all 
fitting was calculated using only statistical errors). 
The inclusive mean pt was extracted for each class as 
(pt) = 2po/{n — 3) (c/. Sec. IVIII for those results). 

In Fig. 121 (left panel) we plot relative fit residuals 
\Jyt Ncvt (data — fit) / V data distributed on yt- The 
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FIG. 2: Left: Relative residuals from power-law fits to pt 
spectra in Fig. The hatched band represents the expected 
statistical errors for STAR data. Right: Exponents n from 
power-law fits to data (solid points) and to corresponding two- 
component fixed-model functions (open circles, see Sec. IVB 
compared to the two-component fixed-model Levy exponent 
12.8 ± 0.15 (hatched band). 



points indicate the actual data positions. The quantity 
plotted insures that the residuals are directly measured 
in units of the r.m.s. statistical error at each yt- These 
relative residuals are then similar to Pearson's correla- 
tion coefficient or relative covariance |0] . Poisson errors 
apply to dn/dyt, whereas the spectra plotted in Fig. ^ 
('data') are of the form \/yt dn/dyt- Thus, a factor ,Jyl 
is required to make the statistical reference uniform on yt 
in these residuals plots. The residuals structure on pt is 
equivalent to that on yt within a Jacobian factor (the fits 
were actually done on pt and the residuals transformed to 
yt for this plot). As noted in the discussion of Fig.llOland 
elsewhere, much of the structure due to hard scattering 
and fragmentation is displaced to small pt in a nonlinear 
way when plotted on pt . 

The large-wavelength residuals in Fig. |21 (left panel) 
exceed the expected statistical error (hatched band) by 
up to 30 X and are similar in form for various rich classes, 
revealing a large systematic disagreement between the 
power-law model and data. The small-wavelength struc- 
ture, mainly attributable to true statistical fluctuations, 
is consistent with expectations (hatched band). The ar- 
gument supporting the power-law model of pt spectra is 
thus shown to fail when tested with high-statistics STAR 
p-p data. 

In Fig. 121 (right panel) we plot best-fit values of power- 
law exponent n vs rich resulting from fits to data (solid 
points) and to the two-component model functions de- 
scribed later in this paper (open circles). The latter 
points and hatched band are discussed in Sec. IXII We 
observe a very strong variation of n with multiplicity. Re- 
duction of n with increasing hard scattering is expected 
in the power-law context, but we find that the physical 
mechanism is different from the theoretical expectation 
(c/. Sec. 1X1. 

We observe very strong disagreement between the 
power-law model function and data, whereas a previous 
UAl (SPPS) analysis reported power-law fits with rea- 



sonable at the same energy The UAl results are 
nevertheless consistent with the present analysis because 
that analysis was inclusive on rich and employed only 20k 
minimum-bias events (?;s 3 x 10^ for the present analy- 
sis). That analysis was therefore statistically insensitive 
to the structures apparent in Fig. [3 Statistics for the 
UAl minimum-bias pt spectrum are comparable to the 
rich = 11.5 multiplicity class in this study, but the latter 
contains about 10 x the hard component in the UAl min- 
bias spectrum. An E735 (FNAL) analysis of spectrome- 
ter data at 0.3, 0.55, 1.0 and 1.8 TeV |li|, including mul- 
tiplicity dependence of spectrum shapes, also obtained 
satisfactory power-law fits to pt spectra. However the ef- 
fective event number was comparable to the UAl study, 
in part because of the reduced angular acceptance of the 
spectrometer relative to the STAR CTB detector, and 
the Pt acceptance [0.15,3] GeV/c was considerably less 
than STAR or UAl, further reducing sensitivity to spec- 
trum shape. Given this exclusion of the power-law model 
we now seek an alternative model which best describes 
Pt spectra from relativistic nuclear collisions. 



IV. RUNNING INTEGRATION 

Running integration provides substantial noise reduc- 
tion for spectrum analysis, thereby improving precision. 
In this section we examine the rich dependence of differ- 
ential and integrated spectra and define alternative nor- 
malization factor ns{fich) and reference function Sq. 



A. Spectrum normalization 

In Fig. 131 (left panel) the spectra from Fig. ^ (right 
panel) are replotted without vertical offsets as spline 
curves for detailed comparison. No assumptions have 
been made about the data, and all spectra integrate to 
unity when extrapolated. The dash-dot curve is reference 
S'o defined in this section. To facilitate the discussion we 
identify three regions on yt separated by the vertical dot- 
ted lines: A = [1.3,1.9], B = [1.9,3.4] and C = [3.4,4.5]. 
The region below yt = 1.3 is outside the pt acceptance. 
Regions A and C are defined such that the curves within 
them are nearly constant relative to one another, whereas 
in region B the differences between curves vary rapidly. 

The trend of the spectra with increasing rich is coun- 
terbalancing changes within A and C: linear decrease in 
A (see inset) and linear increase in C. The relative vari- 
ation in the two regions over the observed rich range is 
quite different: 10% reduction in A and 10 x increase in 
C. Such balancing variations are expected if the yield 
in C increases relative to A with rich, due to the require- 
ment that the normalized spectra must integrate to unity. 
We conclude that with increasing rich additional particle 
yield localized on yt and dominating region C is added 
to the spectrum. 
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The apparent reduction at smaller yt is then a trivial 
effect of the unit-integral condition which can be com- 
pensated by changing the normalization. We normal- 
ize the spectra not by true total multiplicity rich but by 
multiplicity Ug defined such that the normalized spec- 
tra approximately coincide within region A. The vari- 
ation of lower end-point positions with rich is compen- 
sated within errors by normalizing with the linear func- 
tion ns(?ich) = 'ifich (1 — 0.013nc;i) (function ns estimates 
multiplicity n^). The negative term compensates the rel- 
ative yield increase at larger yt. The revised normal- 
ization also facilitates the running integration study de- 
scribed below. 




FIG. 3: Left: Spectra from Fig. (right panel) replotted 
as spline curves and without offsets (solid curves) compared 
to reference (dash-dot curve). Right: Running integrals 
Eq. (0 of extrapolated yt spectra in Fig.^divided by Us/rich 
(solid curves) compared to running integral No{yt) of refer- 
ence So{yt) (dash-dot curve). The ten data curves from bot- 
tom to top correspond to increasing hch G [1, 11.5]. 



Given the results in Fig. |3| (right panel) the natural 
choice for a reference is one which coincides with all data 
curves for yt < 2 and defines a limiting case for the se- 
quence of separated data curves at larger yt- Wc there- 
fore define the reference as the asymptotic limit of the 
yt spectra (or their integrals) as rich — > 0. For reasons 
discussed below we chose as a trial reference the Levy 
distribution 17] 

S'o(mt; /3o, n) = AJ{1 + po (m^ - too)/^)" (2) 

defined on transverse mass rrit and suitably transformed 
to yt- f3o = 1/Tq is an inverse-slope parameter. We 
find that the Levy distribution with optimized parame- 
ters (dash-dot reference curves in Fig. ^ coincides with 
the desired asymptotic form. Determination of param- 
eters n and Pq from the data is discussed in the next 
subsection. Amplitude As{Po,n) is defined by the unit- 
integral normalization requirement on 5*0. 

The running integral of Sq, the dash-dot curve in Fig.O 
(right panel) denoted by iVp, is obtained by replacing 
the curly bracket in Eq. with So{yt), in which case 
n{fich, yt) No{yt) (also, see the legend in Fig.|3|- right 
panel). A^o is thereby defined as the limit as rich 
of the running integrals for the ten multiplicity classes. 
We can obtain a more differential picture by optimizing 
reference curve Sq and subtracting it and its running in- 
tegral Nq from the data. Fig. 0] (left panel) discussed 
in the next subsection reveals the nc^-dependent yield 
increase as a localized structure on yt and is used to op- 
timize 5*0. This differential procedure represents a new 
level of precision in spectrum analysis facilitated by the 
high-statistics STAR p-p data and the running-integral 
technique. 



B. Running integrals and reference So 

To calculate running integrals the measured spectra 
are extrapolated in the pt interval [0,0.2] GeV/c (yt G 
[0, 1.15]) with reference function ris/nch ■ Sq. The ex- 
trapolation is relatively insensitive to the Sq parameters, 
insuring quick convergence of the 5*0 optimization pro- 
cedure described below. The running integral of a yt 
spectrum is defined by 

rvt 

n{hch,yt) = / dy'ty't {l/y'tdn{fich,y't)/dy't} . (1) 
Jo 

In Fig. (right panel) the normalized running integrals 
l/risijich) ■ n{hch,yt) reveal the detailed structure of the 
spectra with much-improved signal-to-noise ratio. We 
observe that the integrals in the right panel indeed nearly 
coincide up to yt ^ 2. Above that point (region B) the 
integrals separate. In region C the integrals all saturate, 
with nearly equal spacings between curves. That result 
provides a first detailed look at the localized (on yt) ad- 
ditional yield which produces the rich dependence of the 
yt spectrum shape. 



C. Optimizing reference So 

In Fig. 0] (left panel) we plot the difference between 
running integrals l/ns{hch) ■ 'ri{hch,yt) of the corrected 
spectra in Fig. (right panel) and reference integral 
^oiut) (the dash-dot curve in that panel). In region B we 
observe a strong localized rich dependence in the yt spec- 
tra. The optimum parameters for Sq are derived as fol- 
lows. Inverse-slope parameter Pq is adjusted to minimize 
residuals in region A of Fig. 21 (left panel). Pq determines 
the average slope of the residuals in that region. Expo- 
nent n then determines the size of the first step in region 
C. n is adjusted so that the first step follows the nearly 
linear trend of rich dependence in that yt interval. Am- 
plitude As{P,n) is determined by the unit-normalization 
requirement for 5*0. 

Changing either Pq or n in Sq does not alter the step- 
wise variation with rich of the data curves in the left panel 
above the first step. That structure is inherent in the 
data and unaffected by the reference choice (c/. Fig. 
right panel, before reference subtraction) . The amplitude 
variation within region C is well represented by Uh /ris = 



yt 



yt 



FIG. 4: Left: Differences between integrals of extrapolated 
yt distributions in Fig. according to Eq. and integral 
No{yt) of soft reference So{yt)- The ten curves correspond to 
rich G [1,11.5] from bottom to top. Right: Distributions in 
the left panel divided by their end-point values at yt = 4.5. 
The dashed curve is the running integral of Ho {cf. Sec. IVt. 
The dash-dot curve is A^o, the running integral of 5*0. 



a rich with a ^ 0.01, where Uh is the coefficient of de- 
fined in the next subsection. That procedure determines 
reference Sq parameters Ag = 20.3 ± 0.1, n = 12.8 ± 0.15 
and To = 0.1445 ± 0.001 GeV |3- 

In Fig. 21 (right panel) the curves are obtained by divid- 
ing the curves in the left panel by their values at upper 
endpoint yt = 4.5 which approximate ratio Uh/ns- Refer- 
ence No{yt) is included in the right panel as the dash-dot 
curve. Comparing A'o to the data integrals it is clear that 
the multiplicity dependence in Fig. ^ cannot be accom- 
modated by adjusting Sq. With the exception of the first 
few rich values (labeled curves) the integrals closely follow 
a common trend: an error function or running integral of 
a gaussian which estimates in a model-independent way 
the running integral of the ric/i-independcnt model func- 
tion iJo(yt) determined differentially in the next section. 



V. DIFFERENTIAL ANALYSIS 

Using running integrals we have defined a precision 
reference for the yt spectra and isolated the Uch depen- 
dence of those spectra relative to the reference. We now 
return to the differential yt spectra and identify an addi- 
tional spectrum component by subtracting the reference 
from the data. The dashed curve in Fig. 0] (right panel) 
(just visible near yt = 2) represents the running integral 
of model function Hq determined in this section. H()[yt) 
models the additional yield at larger yt as a differential yt 
spectrum component. It is already clear from Fig.^that 
the shape of that component is approximately gaussian 
and nearly independent of Uch- 

In Fig. O (left panel) we show the result of subtracting 
reference So{yt) from the yt spectra in Fig. ^ (right panel) 
divided by ng/uch- We obtain the difference distributions 
denoted by Uh/us- H{fich, Vt) (the data points connected 
with dashed curves). Those data represent all rich de- 
pendence of the yt spectra relative to fixed reference 5*0. 



FIG. 5: Left panel: Distributions on yt in Fig. (right panel) 
divided by Us/uch minus reference So{yt)- Dashed curves in 
the left panel and solid curves in the right panel are spline fits 
to guide the eye. The vertical dotted lines enclose region B 
previously defined. Right panel: Distributions H{nch,yt) ob- 
tained by dividing the curves in the left panel by Uh/ris. The 
dashed curve represents hard reference Ho{yt)- The dash-dot 
curve represents soft reference So{yt)- The solid curve under- 
lying the dash-dot curve is an error function !l^ . The hatched 
region estimates the systematic error from subtraction of So. 



The error bars denote statistical errors, applicable also 
to the data in Fig. The two vertical dotted lines en- 
close region B on yt previously defined. H(nch,yt) has 
unit integral by definition, consistent with n^+nh — Uch- 
The shapes of the data curves are well-approximated by 
the unit-integral gaussian reference 



Ho{yt;yt,cryJ = Ah{yt,cryJ ■ exp <^ -- 





yt - yt' 


1 









,(3) 



with Ah = 0.335 ± 0.005, yt = 2.66 ± 0.02 and = 
0.445±0.005. The solid curves represent nh/ng-HQ, with 
best-fit amplitudes nh{fich) /nsifich) plotted in Fig. [7| 
(right panel, sohd dots). Uh is the multiphcity of the new 
spectrum component. The data are generally well de- 
scribed by the model, except for the excursions at smaller 
yt for the smaller rich values. 

Dividing the data in Fig. (left panel) by the cor- 
responding best-fit gaussian amplitudes Uh/us reveals 
the normalized data distributions H{nch,yt) in the right 
panel. Reference So{yt), shown as the dash-dot curve in 
the right panel, is approximately an error function |l9l |. 
The hatched region estimates the systematic error from 
the 5*0 subtraction. Deviations from the Hq model func- 
tion (dashed curve) in that panel represent all the resid- 
ual Hch dependence of the yt spectra, i.e., all deviations 
from the two-component model in Eq. Q below. Those 
deviations are plotted in Fig. (left panel) and discussed 
further in the next section. 

The QCD-based power-law trend p^" expected for 
hard parton scattering would appear in this plotting for- 
mat as a straight line with negative slope equal to the 
exponent or 'power' —n since yt ~ ln(2pt/mo) at 
large pt makes the plot effectively a log-log plot. Out to 
yt — 4.5 or = 6 GeV/c we observe no linear tangential 
departure from gaussian model Hq (dashed parabola) in 
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data H{hch,yt)- 

VI. TWO-COMPONENT MODEL 

The two-component model 0, |23| states that the 
minimum-bias frequency distribution on event multiplic- 
ity from relativistic p-p collisions can be resolved into two 
components, each approximated by a negative binomial 
distribution (NBD) with its own mean and k parameter. 
The two components correspond to events with (hard) 
and without (soft) significant hard parton scatters. That 
concept can be extended to the possibility that the inclu- 
sive pt spectrum shape for hard events is different from 
that for soft events |23| — that the former contains an 
additional spectrum component which we designate the 
hard component, the complement being then the soft 
spectrum component. In that interpretation spectra from 
different multiplicity classes should contain different ad- 
mixtures of the two spectrum components, and the mul- 
tiplicity dependence of the spectrum shape may therefore 
provide a means to isolate those components. 

In this section we examine the two-component model 
in detail. We consider the factorization structure of the 
model function that has emerged from data analysis, we 
examine the residuals structure compared to statistical 
errors and then test the necessity of the fixed-parameter 
model function by fitting the data with all model param- 
eters freely varying. We finally relate all multiplicities in 
the model and show that they form a consistent system. 

A. Two-component model function 

We have analyzed the multiplicity dependence of yt 
spectra from p-p collisions without an a priori model and 
have observed a strong rich dependence whose functional 
forms we now summarize. The two-component model of 
yt spectrum structure can be generally represented by 
the first line of 

1/ytdn/dyt = s{nch,yt) + h{nch,yt) (4) 
= risifich) So{yt) + rihifich) Ho{yt) + . . . , 

with unspecified soft and hard spectrum components 
s{nch,yt) and hifich^yt)- What we have inferred from 
the rich dependence of the measured yt spectra is the 
second line, which represents a factorization hypothe- 
sis with spectrum components modeled by unit-normal 
functions So{yt) and iJo(j/t) independent of rich, ratio 
nhifich)/ risifich) = afich, and constraint rig + rih = rich- 
We suggest that the algebraic model in the second line 
corresponds to the two-component physical model de- 
scribed above and represented by the first line. In the rest 
of this section we consider the quality and details of the 
parameterized model in Eq. and test its uniqueness 
by performing a free fit of the unconstrained model 
functions to the data. 



In the power-law context there is no a priori hypothesis 
for rich dependence: each of the ten multiplicity classes in 
this analysis can be fitted independently with the three- 
parameter model to produce 30 fit parameters. The cor- 
responding residuals are shown in Fig. |2] (left panel) . For 
the two-component model we could in principle have six 
free parameters for each rich, producing 60 fit parameters. 
However, the algebraic model of Eq. |0J (second line) con- 
tains constraints motivated by the requirement of model 
simplicity which greatly reduce the number of indepen- 
dent parameters. 1) The shapes of unit-integral functions 
Soiyt) and Ho^yt) are independent of multiplicity: each 
function is determined by only two parameters fixed for 
all rich- 2) The relative normalization of the two com- 
ponents is nearly linearly proportional to the observed 
multiplicity, as defined by fifth parameter a. Thus, only 
five parameters represent all the data in that model. As 
with the power-law model we compare data to model on 
the basis of relative fit residuals on yt, which provide a 
more differential and direct assessment of fit quality than 
the statistic. 



B. Five-parameter fixed model 

The residuals in Fig. (left panel) correspond to the 
function in the second line of Eq. |0J with five optimized 
parameters held fixed for all rich- Above yt = 2.7 the 
residuals are consistent with statistical fiuctuations ex- 
cept for a few sharp structures with amplitude several 
times the statistical error. Those structures arise from 
the comparatively low statistics of the Monte Carlo sim- 
ulations used for background corrections. The Monte 
Carlo statistical fiuctuations appear in these residuals as 
small- wavelength systematic deviations. 

The prominent residuals in yt < 2.7 for rich = 1-4 (a 
'third component') could represent nontrivial rich depen- 
dence of the soft or hard component or some additional 
physical mechanism. The endpoint values at yt = 4.5 in 
Fig. ^ (left panel) vary linearly with rich to a few per- 
cent (open symbols in Fig. [7| - right panel) , despite the 
substantial nonlinear excursions at small yt of the dis- 
tributions in Figs. 13 and That apparent contradic- 
tion suggests that the prominent residuals may repre- 
sent a change of the hard component at small Uch which 
preserves the linear trend of the integrals. These two- 
component residuals from the five-parameter fixed model 
are otherwise much smaller than the systematic devia- 
tions of the power-law model in Fig. |21 (left panel) with 
its 30-parameter fit, especially in the large- ?/t region 
where the power-law model should be most applicable. 

C. Two-component free fits 

To determine whether the algebraic model of Eq. Q is 
necessary (required by the data), not simply an accident 
of data manipulation, spectra for rich € [1, 11.5] were fit- 
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FIG. 6: Left: Relative residuals between data yt spectra in 
Fig. (right panel) and the two-component fixed parameter- 
ization, for all multiplicity classes. Right; Relative residuals 
between data yt spectra and two-component free fits, also 
for all multiplicity classes. The hatched bands represent the 
expected statistical errors for STAR data. 



ted with the six-parameter function in Eq. Q using 
minimization. Spectra 1/ytdn/dyt were first normalized 
by multiplicity estimator ns{hch) from the fixed param- 
eterization. The coefficients of 5*0 and Hq in the fitting 
function are then ns/fig and nh/fis- The six parameters 
(rig, /3oi nh,yt, o'yt) were freely varied for each rich- 

The residuals from the free fits are shown in Fig. El 
(right panel). The fit residuals are comparable to the 
corresponding fixed-model residuals in Fig. |H| (left panel), 
even though the free fits include six independent param- 
eters for each of ten rich classes for a total of 60 parame- 
ters, compared to the fixed model with only five param- 
eters to describe all ten rich classes. The residuals for 
the smaller hch values show that the free fit attempts to 
minimize the small-j/t structure ('third component') in 
the left panel at the expense of increased intermediate- 
Ut residuals. The effect on the fit parameters is however 
modest, as illustrated in Table |l| 

Table ^ compares the fixed-model parameter values 
(fixed) to the results of the six-parameter free fits (fit- 
ted) for ten rich classes. If the hard-component gaussian 
on ut were not necessary we would expect the fit to 
converge to the soft-component Levy distribution as a 
proxy for the power-law function. The results in Tabled 
indicate that most of the free-fit Sq and Hq shape param- 
eters remain nearly constant within errors across the full 
fich interval. The hard-component gaussian amplitudes 
are definitely nonzero and monotonically increasing, con- 
sistent with the trends in Fig. |S1 (left panel) obtained by 
subtracting So{yt) from the normalized spectra in Fig. ^ 
(right panel). 

Fig. [71 (left panel) shows trends for the two fit pa- 
rameters n and nh/fig which best illustrate the trade- 
off between soft / power-law and hard components of the 
model and the necessity of the two-component model. 
Best-fit values are presented for all rich classes as the 
solid symbols (open symbols are discussed below). There 
are significant systematic deviations of exponent n from 
the fixed-model value (hatched band) which are however 



fitted 


soft component 


hard component 






ris/ris 


ro(GeV) 


n 


nh/fis 


yt 




A. 1 


1 


0.995 


0.145 


11.97 


0.000 






73.3 


2 


1.001 


0.145 


11.78 


0.002 


2.75 


0.500 


40.4 


3 


1.001 


0.145 


11.83 


0.013 


2.75 


0.421 


15.0 


4 


0.996 


0.145 


11.74 


0.025 


2.75 


0.400 


7.36 





0.994 


0.145 


12.60 


0.049 


2.65 


0.427 




6 


1.001 


0.144 


15.63 


0.089 


2.57 


0.450 


1.09 


7 


0.999 


0.144 


15.42 


0.097 


2.57 


0.451 


0.62 


8 


1.005 


0.144 


16.73 


0.115 


2.56 


0.454 


1.18 


9.5 


1.011 


0.143 


16.66 


0.130 


2.56 


0.456 


0.52 


11.5 


0.995 


0.145 


15.69 


0.128 


2.58 


0.460 


1.20 


fixed 


1.000 


0.1445 


12.8 




2.66 


0.445 




error 


0.005 


0.001 


0.15 




0.02 


0.005 





TABLE I: Two-component fit parameters. The line la- 
beled 'fixed' contains the two-component fixed-model param- 
eters. Each fit has v = 2% degrees of freedom. The error 
row applies only to the fixed parameterization. The fit errors 
are generally smaller than those errors for the last five free-fit 
rows. Significant systematic effects are discussed in the text. 
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FIG. 7: Left: Parameters n and nh/fis vs rich from free fits 
(solid symbols) in Table I. The open symbols represent similar 
free fits but with yt held fixed at 2.65 (see text for discussion). 
The bands represent the corresponding two-component fixed 
parameterizations with their stated errors also given in Ta- 
ble I. Right: Ratio nh{hch) /naifich) derived from integrals 
in Fig|l] (left panel) (open squares) and from gaussian am- 
plitudes in Fig. 13 (left panel) (solid dots). The dashed and 
dotted lines have slopes 0.105 and 0.095 respectively. The 
solid curve is described in the text. 



qualitatively different from the trends in Fig. [3 (right 
panel) . The hard-component amplitude rih /fis also devi- 
ates from the linear fixed-model trend, but the trend of 
monotonic increase is even stronger. The hard compo- 
nent appears to be more favored by the free fit than by 
the fixed parameterization. We discuss the systematic 
differences between fixed model and free fit in the fol- 
lowing paragraphs. However, this fitting exercise does 
demonstrate that for almost all rich a two-component 
model is indeed necessary to describe RHIC p-p data. 

The systematic deviations between free fits and fixed 
model in Fig. [7| (left panel) are easily understood. We 
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separately consider rich < 5 and rich > 5 (separated 
by the dotted lines). Generally, there is a strong posi- 
tive correlation between soft-component exponent n and 
hard-component relative amplitude rih/ris originating 
from the requirement to describe the large-yt yield. If 
n decreases the Levy distribution tail rises and the am- 
plitude of the hard-component gaussian amplitude must 
decrease as well to compensate at large yt , and conversely. 
The systematic deviations relative to the fixed model for 
rich < 5 respond to the presence of the 'third compo- 
nent,' which is not a part of the two-component model. 
To compensate for the additional component in the data 
the hard-component amplitude is suppressed and n is re- 
duced by about 10% to provide additional yield from So 
at small yt- The consequence is negative residuals near 
yt = 2.6 in Fig. El (right panel). 

For rich > 5 a different issue arises. In Fig. [51 
(left panel) we have noted previously that the hard- 
component data peaks are skewed (fall off more rapidly 
on the low-yt side) whereas the hard-component model 
function is a symmetric gaussian. The difference is most 
apparent in the running integrals of Fig. 01 (right panel): 
the dashed model curve lies above the data near ~ 2. 
In Fig. (right panel) the hatched region illustrates the 
region of maximum influence of the 5*0 subtraction on 
the hard component. Because the hard-component data 
peaks are asymmetric the 5*0 subtraction at larger yt 
must be reduced by increasing exponent n (the small- 
yt Sq contribution must remain constant to describe the 
spectra there). This requires a compensating increase in 
the hard-component amplitude to fit the larger-j/t part of 
the spectra, and the gaussian model function must shift 
down on yt (by ~ 0.1 or 5 sigma) and the width increase 
slightly (0.01 or 2 sigma) to accommodate the apparent 
increased symmetry of the data hard component. 

To test that description the free fits were redone with 
the gaussian centroid fixed at yt = 2.65. The open sym- 
bols in Fig. d (left panel) show the result. The best-fit 
parameters are now within the error bands of the fixed 
model, with only modest increase in x^/j^ (1.69, 1.07, 
1.43, 0.95, 1.18 respectively for hch — 6, - ■ ■ , 11.5 com- 
pared to the corresponding values in Table P) . The fit 
residuals in Fig. (right panel) appear identical for the 
two cases. We emphasize that the mode (most probable 
point) of the data hard-component peak is near yt = 2.65. 
The downward shift of the model peak in the free fit is a 
consequence of the skewness in the data hard component 
not described by the fixed model but consistent with mea- 
sured fragmentation functions from reconstructed jets. 



D. Two-component multiplicities 

In Sec. IIV Al we adopted a normalization strategy 
which brought all spectra into coincidence in region A 
of Fig. |31 (left panel) by defining multiplicity rig cx rich 
except for a small deviation linear in hch- We then de- 
fined reference function 5*0 as a limiting case of the spec- 



trum rich dependence and isolated a second component 
Hq by subtracting the fixed reference from all spectra. 
The amplitude of Hq relative to the reference is defined 
by ratio rih/ris oc hch- The representation to that point is 
(physics) model independent, derived only from the ob- 
served spectrum rich dependence: the reference is cx hch 
and the second component is cx fi^^. That difference is 
the underlying basis for distinguishing the two compo- 
nents. 

In this section we have identified the two algebraic 
spectrum components with the components of a physical 
model of soft and hard parton scattering and subsequent 
fragmentation to detected particles. We distinguish four 
event multiplicities: 1) the observed multiplicity hch or 
uncorrected number of particles with pt > 0.2 GeV/c in 
the STAR angular acceptance which serves as an event- 
class index, 2) the corrected and pt-extrapolated multi- 
plicity Hch, 3) the 'soft-component' multiplicity Hs and 4) 
the 'hard component' multiplicity Uh, with ris+Uh = Uch- 
We now examine the self-consistency of the multiplicities 
in our two-component model in the context of real spec- 
trum properties, including efficiencies and acceptances. 

Soft multiplicity Us is estimated by function hs{hch) — 
[2.0 ± 0.02(rel) ± 0.2(abs)] hch [1 - (0.013 ± 0.0005) hch]- 
The 1% error apphes to the relative or spectrum-to- 
spectrum normalization relevant to this differential anal- 
ysis, whereas the 10% error applies to the common nor- 
malization of all spectra. As noted, coefficient 0.013 is de- 
termined by requiring that corrected spectra normalized 
by hg approximately coincide within region A of Fig. (21 
(left panel) for all hch- The factor 2 is determined by 
requiring that after correction, extrapolation with 5*0 to 
Pt — and normalization with Uch all spectra in Fig. ^ 
integrate to unity. In the first column of Table ^ devi- 
ations of ris/hs from unity are consistent with the 1% 
error estimate. 

The hard fraction Uh /n^ — a hch is estimated by 
two methods. In the first method we determine the 
gaussian amplitudes required to fit the data distribu- 
tions in Fig. (left panel). Those amplitudes give the 
solid gaussian curves compared to data in that plot 
and are plotted as the solid points in Fig. [3 (right 
panel). The linear trend (dashed line) corresponds to 
slope a — 0.0105 ± 0.0005. The solid curve pass- 
ing precisely through the points is rih{hch)/ns{hch) — 
{(0.0105 r\,,)-io -h (0.005 ni;f)-i°}-i/i°, the errors on 
the coefficients being ±0.0005. The nonlinearity of that 
curve is related to the non-gaussian small-yt structure for 
small values of hch (third component). 

In the second method we note that the distributions 
in the left panel of Fig. 01 are running integrals of data 
distributions in Fig. [3 (left panel). The amplitudes 
of those integrals at end-point yt = 4.5, plotted as 
open squares in Fig. [71 (right panel), also estimate ra- 
tio Uh/ris. They vary nearly linearly (dotted line) with 
slope a = 0.0095±0.0005. Reduction of a from 0.0105 for 
the gaussian amplitudes to 0.0095 for the integral end- 
points results from small deviations of the data peaks 
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from the Hq gaussian model at small yt evident in Fig.O 
The data are slightly skewed in a manner consistent with 
measured fragmentation functions. The model gaussians 
are matched to the data at and above the data peak mode 
or most probable point. The integral of any data peak is 
therefore expected to be slightly less than that of the cor- 
responding model function. Both methods suggest satu- 
ration of the hard-component amplitude at larger ric/i. 

Consistency of the soft and hard multiplicity estima- 
tors within the two-component model can be established 
by the following argument: Tracking inefficiencies pro- 
duce the same fractional changes for all rich and are repre- 
sented by factors and en for soft- and hard-component 
yields. The corrected spectra are extrapolated to pt — 
with soft model 5*0. The fraction of So falling above 
Pt = 0.2 GeV/c (within the pt acceptance) is represented 
by 7. The hard component identified in this analysis falls 
entirely within the pt acceptance. The observed multi- 
plicity is then given by rich = l^s^s + ^h^h, whereas 
the corrected and extrapolated spectra integrate to true 
multiplicity rich = rig+rih- The expression for rich above 
can be rearranged to solve for rig in the first line below, 

Us — — ^ < 1 — ■ a hch (■ predicted (5) 

hs = 2nc?i {1 — 0.013 fic/i} observed 

whereas the second line is the estimator inferred from 
the data. By integrating reference So we determine that 
7 = 0.7: 70% of the reference spectrum is within the 
acceptance pt > 0.2 GeV/c. Tracking efficiencies eg and 
eh are both approximately 70%, and we have determined 
from the data (running integrals) that a ^ 0.0095. We 
therefore have 1/76^ ~ 2 and eh/'yeg ■ a ~ 0.0135, es- 
tablishing the consistency (predicted ^ observed) of the 
two-component multiplicities. Coefficient 0.013 is iden- 
tified as a/7, and the trend of is defined by ratio 
nh/ris = enrich- We thus close the circle, demonstrat- 
ing quantitatively how increase with rich of the hard- 
component contribution to the spectrum forces to de- 
crease relative to hch in compensation, why hs contains 
the negative term and what its magnitude must be. 

VII. (pt) SYSTEMATICS 

Another aspect of the two-component model is the 
variation of ipt) (inclusive meanpt) with hch- Estimation 
of {pt) for spectra with incomplete pt acceptance requires 
either a model fit or direct integration of data with ex- 
trapolation. The power-law function for pt distributions 
l/ptdn/dpt = A/{l+pt/poT, with lj>t) = 2po/{n ~ 3), 
has been used previously to extract (pt) values from 
corrected pt spectra {pt) can also be determined 

by direct integration of the experimental pt spectra, 
with extrapolation to pt = by a suitable model func- 
tion. Finally, the two-component fixed-model function 
obtained in this analysis can provide a parameterization 
of {pt){nch)- 



The running multiplicity integral n{hch,yt) is defined 
by Eq. with the data extrapolated over pt G [0,0.2] 
GeV/c by reference ns{nch) So{pt) - Running integral 
Pt{hch,yt) can also be defined for transverse momentum 
Pt by including an extra factor Pt{yt) iu the integrand of 
Eq. IU. The ratio {pt){hch,yt) = Pt{hch,yt)/n{nch,yt) is 
then a function of yt for each value of hch, and (pt){hch) 
is the limit of that function as yt ^ 00. {pt)(hch) is thus 
determined by direct integration of pt or yt spectra. 

A changing mixture of soft and hard components may 
cause (pt) to vary with rich- The (pt) values for individual 
components arc obtained by direct integration of model 
functions 5*0 and Ho'- {pt)soft = 0.385 ± 0.02 GeV/c and 
{ptihard — 1.18±0.01 GcV/c. A two-compoueut analytic 
expression for {pt) is then given by 

lj>,)[fich) = (o.385^^^^ + i.is^^^^l GeV/c, (6) 
L nch nch J 

with Uh/us = a hch and Us + Uh = Uch- 
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FIG. 8: {pt){hch) derived from the two-component Ho gaus- 
sian amplitudes (solid dots), from the running integrals (open 
triangles) and from power-law fits to STAR and UAl data 
(open circles, triangles). The solid and dotted lines corre- 
spond to Eq. Q with a — 0.0095 and 0.015 respectively. 

In Fig. IHl (pt) (rich) values inferred from power-law fits 
to corrected STAR spectra are represented by open cir- 
cles, consistent with a 200 GeV UAl power-law analysis 
plotted as solid triangles 8], but inconsistent at smaller 
hch with the two-component result from this analysis 
plotted as solid points. The two-component data were 
obtained with the Uh/us values plotted as solid dots in 
Fig. 13 (right panel). The solid line represents the two- 
component analytic expression for (pt) in Eq. © with 
a = 0.0095. The UAl results for 900 GeV % are plotted 
as open triangles. The dotted line corresponds to Eq. ® 
with a =^ 0.015. (pt) values obtained by direct integration 
of the extrapolated spectra are represented by the open 
squares. The hatched region represents the common un- 
certainty in all means due to uncertainty in the particle 
yield in pt < 0.2 GeV/c. 

The {pt){hch) values in Fig. |S1 obtained by direct in- 
tegration of extrapolated spectra provide the best esti- 
mate of the physical trend. The results at 200 GeV for 
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direct integration, the two-component model and power- 
law fits are consistent within errors for rich > 4. The 
notable deviation of the power-law results from the two- 
component linear trend for hch < 5 can be explained 
by the third-component structures at small yt and small 
■hch in Fig. (left panel). Those structures strongly influ- 
ence (bias) extrapolation of the power-law function into 
the unmeasured region in pt < 0.2 GeV/c so as to over- 
estimate the inferred yield there (nominally 30% of the 
total spectrum). The overestimate at small produces a 
sharp reduction of ipt){fi-ch) values inferred from power- 
law fits. The additional yield at small yt in Fig. [3 itself 
corresponds to {pt) ^ 0.4 GeV/c, and thus cannot phys- 
ically lower the composite {pt) below {pt)soft = 0.385 
GeV/c. These {pt) results demonstrate that the UAl 
data are sensitive to the small-j/f and small-nc/i struc- 
tures revealed in this analysis when the more integral 
spectrum measure (pt) is used. 



VIII. ERRORS 

The statistical errors for the basic yt spectra in Fig. 
are best indicated by the error bars on the difference 
distributions of Fig. |31 (left panel). That figure also 
compares the point-wise statistical errors to the hard- 
component structure inferred in this study, which is sta- 
tistically well determined for all rich classes. Monte Carlo 
calculations of background corrections with full detector 
response simulation are computer intensive. Because of 
limited statistics the statistical fluctuations in the Monte 
Carlo data used for background corrections are injected 
into the corrected data spectra as visible systematic er- 
rors: long-wavelength systematic error is reduced at the 
expense of increased short-wavelength random 'system- 
atic' error. Those errors are apparent as the nonstatis- 
tical short-wavelength structures in Figs. [3 and El The 
systematic uncertainties in the corrected spectra can be 
divided into ric/i-dependent and rich-independent uncer- 
tainties. 

nch-independent systematic uncertainties include un- 
certainties in the corrections for tracking efficiency, back- 
grounds (mainly weak decays) and momentum resolu- 
tion. Systematic spectrum corrections for this analysis 
were 20% or less, except for the lowest two pt bins where 
they increased to 40%. Statistical errors for the system- 
atic corrections were typically less than 1% (except as 
noted above for the background corrections). We esti- 
mate the uncertainties in the systematic corrections as 
10% of the correction values. The total uncertainty for 
the systematic corrections is then less than 2% above pt = 
0.4 GeV/c. The UAl corrected ric/i-inclusive spectrum 
for 200 GeV p-p collisions agrees with the correspond- 
ing inclusive spectrum from the present analysis at the 
2% level. 

nch-dependent systematic errors could result from rich- 
dependent tracking inefficiencies. However, track detec- 
tion and Pt measurement in this analysis required no 



reference to other tracks or a fitted event vertex, thus 
minimizing any rich dependencies. In effect, each track 
was treated in isolation independent of its relationship 
to any event, except for the timing requirement with the 
CTB. The tracking efficiencies for low- multiplicity (1-4) 
and high-multiplicity (> 4) events integrated over the pt 
acceptance were found to be consistent to 3%, with a 
1% statistical error. We take that as an estimate of the 
ric/i-dependent systematic uncertainty. 

The main source of systematic uncertainty in the shape 
of the hard-component structures isolated in Fig.[Slis the 
definition of Sq as the lowest element of the regular se- 
quence in Fig. ^ (left panel). So is a rapidly-decreasing 
function in the interval yt = 1.6-3. The main effect of 
varying either /3o or n in 5*0 is to change the magnitude of 
Sq in that interval, shape changes being secondary. It is 
consistent within the two-component context to require 
that 1) component H{nch,yt) be non-negative, placing 
an upper bound on 5*0 in Fig. |S1 and 2) that any rich- 
independent aspect of the distributions in Fig. |Slbc min- 
imized, determining a lower bound. Those criteria place 
stringent constraints on 5*0 already in yt ~ 1.6-2, limit- 
ing systematic offsets at yt = 2 to ±0.002, the allowed 
range rapidly decreasing above that point according to 
the 5*0 curve in Fig. ^| (right panel). The systematic un- 
certainty estimate corresponding to those trends is rep- 
resented by the hatched region in Fig. (right panel). 

The nonstatistical power-law fit residuals in Fig.|21are 
as much as thirty times the statistical error. One of 
the findings of this study is that the power-law model 
function is inappropriate for these pt spectra. System- 
atic uncertainties for the fit parameters are therefore not 
meaningful. 

The fitting uncertainties for the fixed-model parame- 
ters are given at the bottom of Table Those uncer- 
tainties are meaningful relative to the fitting procedure 
defined in the two-component model context. The ability 
of that model to describe the data is apparent in Fig. El 
(left panel) . The only significant residuals correspond to 
a low-yt spectrum element (for hch = 1-4) deliberately 
omitted from the two-component model. One source of 
systematic uncertainty in those parameters is whether 
the fixed-model prescription forces a certain result by ex- 
cluding some other which may better describe the data. 

To test that possibility a free fit with all model 
parameters varying was conducted. The difference in 
the two cases is summarized in Table U and Fig. |7| (left 
panel). In particular, there are substantial differences 
in the Levy exponent and the hard-component ampli- 
tude for the free fit depending on whether the position 
of the hard-component gaussian is constrained or not. 
When the gaussian position is constrained the free fit and 
the fixed model agree within the systematic uncertain- 
ties in the latter. The differences in the unconstrained 
fit are traced to significant departures of the shape of the 
hard component data peak from the symmetric gaussian 
peak shape: the model function could be further refined 
by adding a skewness (expected for fragmentation func- 
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tions) to improve the stability of the fits. However, it 
is not our purpose to develop a complex representation 
of yt spectra, but rather to demonstrate the essential 
two-component aspects of the spectra with the simplest 
possible model function. The differences in fit param- 
eters in Fig. (left panel) can therefore be taken as a 
generous estimate of the systematic uncertainty in the 
fixed-component parameterization. 



IX. IDENTIFIED PARTICLES 

Model functions 5*0 and Hq derived from this analy- 
sis of unidentified particles represent physical spectrum 
components S and H for several hadron types, mainly tt, 
K and p. Two questions emerge: 1) to what extent do 
and Hq correspond to individual hadron types, and 2) to 
what extent does the rich dependence of the pt spectrum 
truly separate two physical components S and H7 The 
soft component of one hadron species may have signif- 
icant rich dependence which could be misinterpreted as 
the hard component of another species, or of the combi- 
nation of unidentified hadrons in this study. 

We can obtain some answers to those questions from 
ric/i-inclusive spectrum studies of identified hadrons. pt 
spectra for -/F/vtv = 200 GeV p-p collisions have been 
measured for identified pions, kaons and protons '2^. 
Because rich ^ 1 and the pt acceptance was [0.3,3] GeV/c 
for that analysis the measured multiplicity-inclusive pt 
spectra are reasonably described by Levy distribution 5*0, 
especially the kaon and proton spectra. The common 
Levy exponent for the three species is n = 16.8 ± 0.05, 
compared to n = 12.8±0.15 measured in this analysis for 
unidentified hadrons. The slope parameter for identified 
pions is T = 0.145±.001 GeV, whereas for both kaons and 
protons T = 0.23±0.005 GeV, compared to T = 0.1445± 
0.001 GeV for unidentified hadrons in this analysis. 

The trend of 5*0 with hadron species is easily under- 
stood. Addition of the 'hotter' K and p spectra to the 
'cooler' pion spectrum flattens the unidentified hadron 
composite at larger pt, reducing the exponent of 5*0 to 
n = 12.8. At smaller pt the pion fraction dominates the 
composite spectrum, and the unidentified-hadron slope 
parameter is the same as the pion slope parameter. The 
effect of the heavier hadrons on the composite spectrum 
is mainly to reduce the Levy exponent from the larger 
physical value common to all three hadron species. 

Information on the rich dependence of pt spectra for 
identified particles is limited. A preliminary analysis of 
Kq and A pt spectra up to 4 GeV/c ||23(| suggests that 
the rich dependence of both spectra can be described by 
a modest (5%) reduction of n with increasing rich- That 
trend can be compared to the free fit results for Sq 
in Tabled as shown in Fig. [7| (left panel): n increases by 
about 25% over the measured rich range. That increase 
is traced to an attempt by the model to accommodate a 
skewness of the hard component in the data, not a true 
variation in the soft component. 



X. PYTHIA MONTE CARLO 



A similar analysis of p-p collisions from the Pythia 
Monte Carlo [2j reveals substantial deviations from 
data. We studied default Pythia- V6. 222 and Pythia 
'tune A' (increased initial-state radiation and multiple 
soft parton scatters relative to the default) with parame- 
ters derived from studies of the underlying event in trig- 
gered jet events In Fig. El (left panels) we show 
Pythia pt spectra normalized to unit integral and soft 
reference 5*0 (dash-dot curves) determined by the same 
criteria applied to STAR data. Those plots can be com- 
pared to Fig. O (left panel). In Fig. |51 (right panels) 
we show the results of subtracting reference 5*0 from the 
normalized spectra in the left panel divided by rig /rich- 
Those plots can be compared to Fig. [S] (left panel). The 
dashed curves are the hard component Hq for STAR data 
divided by 10 to provide a reference. 




FIG. 9; Two-component analysis applied to Pythia Monte 
Carlo data with the same multiplicity classes as for STAR 
data. The left and right panels may be compared with Figs. 
1^ (left panel) and |^ (left panel) respectively. Dashed curves 
in the right panels represent the STAR data hard component 
for rich ~ 11 (Ho/ 10). Dash-dot curves in the left panels 
represent soft component optimized for each Monte Carlo 
configuration: Pythia V6.222 default with parameters To = 
0.147 GeV and n = 23 (upper panel); Pythia Tune A with 
parameters To = 0.137 GeV and n = 14 (lower panel). 

The So parameters for Pythia- V6. 222 in the upper 
panels are To = 0.147 GeV and n = 23. The large 
value of n implies that the Pythia soft component is 
nearly Maxwell-Boltzmann, in sharp contrast to RHIC 
data. The exponent is strictly limited to a large value 
by the Pythia data for yt < 2.5. The Sq parameters for 
Pythia tune A in the lower panels are Tq = 0.137 GeV 
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and n = 14. The smaller value of n is comparable to the 
value 12.8 observed for RHIC data. 

The hard-component yield for Pythia is generally a 
factor of two to three less than the data (most apparent 
above yt — 2.7), broader and peaked at a smaller value 
of yt- Pythia- V6. 222 shows a saturating of the hard- 
component amplitude with increasing hch, whereas Tune 
A shows a more uniform and significantly greater rate 
of increase. The large gaussian-shaped offset common to 
all curves and centered at yt '-^ 2 is also not observed 
in the data. That structure cannot be accommodated 
by the Levy distribution. The two Pythia Monte Car- 
los thus exhibit some features which agree qualitatively 
with experimental data but are quantitatively different. 
Tune A is closer to data than the default for soft and 
hard components, but the ric/i-independent gaussian off- 
set near yt = 2 persists and is not observed in the data. 



XI. DISCUSSION 

A description of p-p collisions in terms of soft and hard 
components is natural at RHIC energies where signifi- 
cant hard parton scattering occurs but the underlying 
event (2^] is still relatively simple. The two-component 
model of nuclear collisions can be applied to 1) the event- 
frequency distribution on rirh (two or more negative- 
binomial distributions) [lOillflj 2) the dependence of (pt) 
on rich lEQil^i 3) triggered jet correlations on (77, 0) 
(correlations from soft and hard event classes) and 4) 
the rich dependence of the pt or yt spectrum shape . 
The common theme is the relation of hard parton scat- 
tering to event multiplicity in the context of a 'soft' un- 
derlying event. This paper emphasizes analysis type 4) 
- study of the rich dependence of the spectrum shape on 
transverse momentum pt and transverse rapidity yt- 

Model functions So{yt) and _ffo(2/t) in Eq- Q) can be 
viewed as the lowest-order elements of a pcrturbative ex- 
pansion of the spectrum shape. Multiplicities ns(rich) 
and rihijich) can be interpreted as estimating the mean 
numbers of soft- and hard-component particles per event 
for a given rich- The claim of simplicity for the two- 
component fixed model is supported by the small number 
of parameters, the simplicity of the model functions, the 
demonstration of necessity in Sec. IVI CI and the demon- 
stration with residuals plots that there is no additional 
information in the spectra (aside from the small-j/f 'third 
component' which may represent additional physics). 

We cannot rule out additional components or changes 
in the shapes of physical components S and H. Each 
should be rich-dependent at some level, but the present 
analysis indicates that within the observed rich interval 
any such dependence is near the level of statistical error. 
A change in S is suggested by the n dependence of the free 
fits in Fig. (left panel). However, that behavior may 
simply be due to a coupling of soft and hard amplitudes 
in the free fit, with no physical significance. 

A significant change in H is expected at larger rich 



based on known jet physics: larger fragment multiplici- 
ties are produced by more energetic partons, with frag- 
ment distributions shifted to larger yt Thus, the 
mean and width of H should increase with rich at some 
point, but such changes are not observed beyond statis- 
tics within the yt and rich acceptances of this study. Ap- 
parently, the multiplicity increase in this analysis is dom- 
inated by increased frequency of events with a single hard 
scattering within a multiplicity class rather than bias to- 
ward more energetic partons. That scenario is consistent 
with the two-component model of [lol| . 

The soft-component Levy distribution 5*0 = As/{1 + 
Poirnt — mo)/n)" is similar in form to power-law 
function A/{1 + pt/po)"^. However, the physical in- 
terpretations are quite different. The Levy distribu- 
tion describes a nominally exponential function with a 
control parameter {e.g., slope parameter) which under- 
goes gaussian-random fluctuations. Inverse exponent 1/n 
then measures the relative variance (Tp/^Q of the control 
parameter jS^I . In the limit 1/n ^ the Levy distribu- 
tion on rrit becomes a true Maxwell-Boltzmann distribu- 
tion. Those properties suggested the Levy distribution 
as a reference function for this analysis. Ironically, the 
'power-law' function in the form of a Levy distribution 
describes the soft component, not hard parton scattering. 
The Levy parameters can be interpreted in the context 
of an ensemble of hadron emitters with random trans- 
verse speeds, thermal radiation from moving sources as 
described by the Cooper- Frye formalism [23 ■ The ex- 
pected QCD hard-scattering power-law trend is not evi- 
dent in the data out to pt ^ 6 GeV/c. 

In Fig. 121 (right panel) we plot exponent n values from 
power-law fits to data (solid points) and to the two- 
component fixed model (open circles) for the full range 
of fich- The latter procedure simulates a power-law fit 
to data with no small-j/t excursions or 'third component' 
and illustrates the effect of those features on the expo- 
nent. The range of variation of the power-law exponent, 
in contrast to the two-component fixed model, and the 
substantial effect of the 'third component' further illus- 
trate that the power-law parameterization is sensitive to 
aspects of spectra inconsistent with its theoretical moti- 
vation, making fit results difficult to interpret physically. 

The gaussian shape of Ho^yt) inferred from this anal- 
ysis can be compared with fragmentation functions from 
jet analysis of p-p, e-p and e-e collisions plotted on loga- 
rithmic variable = ln|io,et / P franm.ent 1 1 which also have 
an approximately gaussian shape [2g | explained in a QCD 
context as the interplay of parton splitting or branching 
at larger pt and the nonperturbative cutoff of the branch- 
ing process at smaller pt due to gluon coherence [sol Isif . 
The gaussian parameters are predicted by the pQCD 
MLLA (modified leading-log approximation) 'd2\. The 
hard component obtained in this analysis then represents 
not fragmentation functions from reconstructed large- i^t 
jets but rather the average of a minimum-bias ensemble 
of fragmentation functions dominated by low-Q^ parton 
scatters [Q < 10 GeV). In that context Hq represents 
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minimum-bias partons dominated by minijets |33l| . A 
previous study of small- i?t clusters in 200 GeV p-p colli- 
sions suggested that semi-hard parton scattering or 
gluon radiation from projectile constituent quarks could 
produce substantial small-p* structure in hadron spectra 
similar to the hard component of this study. 

A recent analysis of pt spectra in the interval 0.3 - 
10 GeV/c for identified particles in p-p and d-Au colli- 
sions [23 used the relativistic-rise particle identification 
scheme to extend the spectra with very good statistics 
to large pt- That paper compared the spectra to several 
NLO pQCD calculations and compared the nit spectra of 
pions and protons. It concluded that there is a transition 
region from soft to hard particle-production processes at 
Pt ^ 2 GeV/c in inclusive particle production, which 
would appear to contradict the present results. How- 
ever, the identified-particle spectra in that study below 
Pt = 2.5 GeV/c are from a previous study |22| | in which 
the point-to-point systematic errors and the statistical er- 
rors are quite large, the latter due to the small acceptance 
of the prototype ToF detector. The ToF-based studies of 
multiplicity-averaged p-p collisions are therefore not sen- 
sitive to the hard-component structure reported in this 
paper, the great majority of which falls below 2.5 GeV/c. 
The present study takes a new approach by comparing 
large-statistics inclusive-hadron spectra in several multi- 
plicity bins. Since the hard component is relatively en- 
hanced in high-multiplicity events we are able to extend 
our investigation of the hard component to low pt by 
studying the trend of that enhancement. 

The relative frequency of hard scatters in p-p collisions 
is described by the fifth model parameter a ^ 0.01, repre- 
senting the nearly-linear dependence of rih /rig on rich- We 
relate the hard-component amplitude to the frequency of 
hard collisions (/ = number of hard collisions per NSD 
p-p collision) as nh{nch) = ahchns{nch) = f{nch) ■ "mj, 
with mean true event multiplicity hch = 2.5 in one 
unit of pseudorapidity and mean minijet multiplicity 
n„ij ~ 2.5 ±1 36!j. We then estimate the observed fre- 
quency of hard scatters in %fs = 200 GeV p-p collisions 
as / = nh/fimj — 0.012 ± 0.004 observed hard scatters 
per NSD p-p collision per unit of pseudorapidity. In that 
interpretation multiplicity rich serves as a 'trigger' for 
hard parton scattering, determining the fraction of hard- 
scattering events in a given multiplicity class and thus 
the relative amplitude of the hard spectrum component. 

Model functions and Hq on pt and yt are summa- 
rized in Fig. ^1 which can be compared with Figs. ^ 01 
and El i?o/9 nh/rig = 0.11 is compared to data for 
f^ch = 11.5 and illustrates the role of the hard component 
in the measured spectra with sufficient amplitude to be 
visible in a linear plotting format (right panel). Simi- 
larly, 7?o/140 =J> nil/ Us = 0.007 is compared to data for 
nch = 1- Those coefficients are consistent with the mea- 
sured hard-component gaussian amplitudes for hch = 1 
and 11.5 (c/. Fig.Q- right panel). 

Collisions in the event ensemble containing at least one 
semi-hard parton scatter within the detector acceptance 
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FIG. 10: Decomposition of inclusive pt spectra into a soft 
component represented by a Levy distribution on mt and a 
hard component represented by a gaussian on yt- Dashed 
curves ffo/9 correspond to data for hch = 11.5, while dotted 
curve Ho/IAQ correspond to data for hch = ^ {cf- Fig.0. The 
dash-dot curves are soft reference So, and the solid curves are 
the totals of soft and hard components for the model. The 
dotted curve in the right panel estimates the shape of the 
inclusive yt distribution for those p-p collisions containing at 
least one minimum-bias hard parton scatter (hard events). 



should have similar yields of soft and hard components 
(assuming an average minijet multiplicity of 2.5). The 
average yt spectrum for such hard events is illustrated 
by + Hq, shown as the dotted curve in Fig. ^| (right 
panel). We cannot isolate such hard events in an unbi- 
ased manner, but we can infer their structure by extrap- 
olating the nch trends determined in this analysis. 

The left panel of Fig.^]indicates the loss of visual sen- 
sitivity to spectrum structure when spectra are plotted 
on Pt ■ The hard component can appear to be a continua- 
tion of the soft component, whereas in the right panel the 
two components are clearly separate functional forms, yt 
provides a more balanced presentation of structure result- 
ing from hard-scattered parton fragmentation, yet does 
not compromise study of the soft component, which is 
well-described by a simple error function on yt jl9l] . The 
transverse and longitudinal fragmentation systems un- 
dergo similar physical processes and should therefore be 
compared in equivalent plotting frameworks. Just as j/z 
is preferred to Pz we prefer yt to pt ■ 



XII. SUMMARY 

In conclusion, we have studied the event multiplicity 
nch dependence of high-statistics transverse momentum 
Pt or transverse rapidity yt spectra from p-p collisions at 
= 200 GeV. We have determined that the 'power-law' 
model function fails to describe the spectra for any nch, 
exhibiting large nonstatistical deviations from data. An 
earlier UAl study reporting satisfactory power-law fits 
to data seems contradictory. However, it is statistically 
consistent with the present study because the UAl data 
were derived from a much smaller event sample. We 



16 



have analyzed the shapes of the spectra with a running- 
integral technique and determined that the spectra can 
be described precisely by a simple five-parameter model 
function. The algebraic model can in turn be related to 
a two-component physical model of nuclear collisions. 

The power-law function motivated by pQCD expec- 
tations for hard parton scattering better describes the 
soft component in the form of a Levy distribution on rut 
(two parameters). We observe for the first time that the 
hard component is well described by a gaussian distribu- 
tion on transverse rapidity pt, with shape approximately 
independent of multiplicity (two parameters). The hard- 
component multiplicity fraction increases almost linearly 
with event multiplicity (the fifth model parameter). A 
detailed comparison of (data — model) residuals from the 
two-component fixed model and from free fits with all 
two-component model parameters varied confirms that 
the two-component fixed model is required by the data. 

The hard component may represent fragments from 
transversely scattered partons. The shape is consistent 
with fragmentation functions observed in LEP and PE- 
TRA e+-e~ and FNAL p-p collisions. The stability of 
the hard-component shape with event multiplicity sug- 
gests that a gaussian distribution on yt is a good rep- 
resentation of minimum-bias parton fragments. The 
relative abundance of soft and hard components at any 
yt of course depends on yt and rich, but most of the hard- 
component yield falls below 2.5 GeV/c. There is evidence 
for a small but significant third component at smaller yt 
and smaller rich- Comparison with the Pythia Monte 
Carlo reveals qualitative differences from data. 
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APPENDIX A: SYMBOL DEFINITIONS 

Below is a list of symbols and their definitions as used 
in this paper. 

yt'. transverse rapidity, replaces transverse momen- 
tum pt to provide improved visual access to 
fragment distributions 

hch- observed event multiplicity in the detector ac- 
ceptance, also the event-class index 
n'^f^: efficiency- and acceptance-corrected multiplic- 
ity in the detector acceptance 
rich- corrected and P(-extrapolated or 'true' multi- 
plicity in the detector angular acceptance 
Us', soft-component multiplicity in the acceptance 
fis- particular function of rich used to estimate rig 
rih- hard-component multiplicity in the acceptance: 

ns + nh = rich 
a: hard-component coefficient: nh/ria ^ ahch 
Sq: unit-normal functional form of the soft compo- 
nent (Levy distribution on rrit) 
Nq: running integral of soft reference So 

Hq: unit-normal functional form of the hard compo- 
nent (gaussian distribution on yt) 

A, po, n: power-law model parameters 

As, /?o, n: soft-component Levy distribution param- 
eters, l/Po = Tq, the slope parameter 
-^h, yt, <^yt'- hard-component gaussian parameters 
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